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^ ; Abstract 

^ ■ This paper considers the single factor Heath-Jarrow-Morton model for the interest rate curve 
■ with stochastic volatility. Its natural formulation, described in terms of stochastic differen- 
00 . tial equations, is solved through Monte Carlo simulations, that usually involve rather large 
computation time, inefficient from a practical (financial) perspective. This model turns to be 
0_( \ Markovian in three dimensions and therefore it can be mapped into a 3D partial differential 
O ' equations problem. We propose an optimized numerical method to solve the 3D PDE model 
P5 ■ in both low computation time and reasonable accuracy, a fundamental criterion for practical 
^ . purposes. The spatial and temporal discretization are performed using finite-difference and 
CT! Crank-Nicholson schemes respectively, and the computational efficiency is largely increased 
performing a scale analysis and using Alternating Direction Implicit schemes. Several nu- 
merical considerations such as convergence criteria or computation time are analyzed and 
discussed. 
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1. Introduction 

In quantitative finance, the interest rate curve has been intensely studied and modeled in 
terms of stochastic differential equations (SDE), assuming for this curve a temporal evolution 
which satisfies the non-arbitrage opportunity in a complete and efficient market H]- Such 
restriction implies that the model at hand needs to be calibrated by the market prices of 
the most liquid instruments, bond prices being a representative example. These in turn are 
decomposed in more elementary units, the so called zero coupon bonds. These are, roughly 
speaking, financial instruments that pay one unit of currency at a certain future date (ma- 
turity) . The price of these latter instruments will characterize the average behavior of the 
interest rate curve On the other hand, the fiuctuations of the interest rate curve with 
respect to its average, the so called volatility, is quantified through market instruments such 
as caps and floors jJ] . A cap pays the difference between a certain rate and a certain prespec- 
ified level (strike) if this difference is positive, while the floor pays off the difference between 
the strike and the value of the rate, if positive. The relation that links the cap or floor 
premium with the volatility is the well known Black-Scholes formula 

Interested in the calibration of zero coupon bond prices together with cap/floor prices, 
in this paper we address the Heath-Jarrow-Morton model (HJM) which in its natural for- 
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mulation is described in [3|. The main interest of the HJM approach is set on the fact 
that it provides a broad mathematical formulation [cj, where most of market observed fea- 
tures can be taken into account. This wasn't the case of preceding seminal models such as 
Black-Derman-Toy model or Hull- White models where the model parameters were very 
difficult to calibrate in terms of market observed patterns. Furthermore, those models only 
incorporated one source of stochasticity, and therefore the only possible mode in the interest 
rate curve was the parallel movement. Only lognormal or normal statistical distributions of 
the short rate were possible. These distributions for the short rate arise in distributions for 
different maturities of the libor rates, which are far away from those implied in the markets. 
Conversely, HJM models are not afflicted by these drawbacks. The HJM framework is 
general, in the sense that many previous models describing the term structure of interest 
rates can be understood as particular cases of a HJM model, that in turn can incorporate as 
many risk factors as needed in order to accurately describe the evolution of the rate's curves. 
The formulation of the HJM can be extended to incorporate several stochastic factors; hence, 
the simulated interest rates curve movements could include deformation modes that changes 
the initial slope and convexity of this curve, in order to describe the covariance and auto- 
correlation structures present in the time behaviour of the interest rate curve. Additionally, 
the probability distribution function of interest rates can be exogenously defined by means 
of local volatility functions (of deterministic or stochastic nature) in order to match option 
prices as quoted in the market. 



The major drawback of HJM models is set on the continuous nature of its state variables 
(the continuous time structure of forward rates), what leads to an infinite amount of state 
variables [6|]- Therefore, in general, these models are non-Markovian in finite dimensions, and 
thus can only be solved by using fitted Monte Carlo techniques which are eventually slow and 
computationally delicate. This is of course problematic for practical purposes, since financial 
industry requires models that can be integrated in real time with reasonable accuracy. Quite 
interestingly, HJM models can be in some particular situations transformed into a low order 
Markovian system 0, H, 0, [lol, 11|, and therefore can be subsequently mapped to a partial 
differential system due to the well known Feynman-Kac formula |l2[|. Numerical methods for 



solving partial differential systems can consequently apply |13l. |14||. Although integration of 
PDE is typically faster than Monte Carlo simulations of the associated SDE system, in prac- 
tice one needs additional optimization methods for calibration (conjugate gradient, genetic 
algorithms), thus requiring a large number of evaluations of the model. To find out fast and 
efficient model evaluations (numerical solvers) is therefore demanding in financial industry. 



In this paper, we propose a full numerical methodology to optimize such issues. Amongst 
the plethora of different HJM models, we address here a single factor HJM model with 
stochastic volatility. This choice is justified on practical reasons (this is a realistic enough 
model which is actually used in the financial industry [l5[ [l6j) and can be argued as it fol- 
lows: (i) single factor: Principal Component Analysis of the covariance matrix associated to 
historic data of the interest rate curve suggest that the first eigenvalue has a weight of 85% 



17|, Il8|. (ii) stochastic volatility: The preliminary HJM models used to model the stochastic 



behavior of the interest rate market by a single Wiener process that drives the forward rate 
processes (that is to say, there was a single source of stochasticity in the models). In the last 



2 



years, some authors have introduced [19|, according to empirical evidence, a new source of 
stochasticity in the description of the volatihty evolution, leading to the so called stochastic 
volatility HJM models |2o|- 



Amongst all the possible stochastic volatility functions, we will focus on those which 
are separable, that is, that can be factorized in the product of a stochastic function and a 
time dependent deterministic function and another deterministic function that depends on 
the maturity. Again, this choice isn't random, much on the contrary, the resulting HJM 
model turns to be Markovian in only three dimensions, and can therefore be mapped and 
evaluated within a PDE framework. We will present a complete numerical integration of 
the model based on finite differences schemes and Crank-Nicholson temporal integrators. In 
order to drastically improve the computational efficiency of the method (that is, reaching fast 
computation time while preserving good accuracy), we employ a scale analysis for the mesh 
optimization and, for the first time in financial models (as far as we are aware). Alternating 
Direct Implicit (ADI) schemes, techniques borrowed from computational fluid dynamics. The 
rest of the paper goes as follows: Section H] focuses on the interest rate market, defining the 
zero coupon bond as the trading derivative. Then, section [3] presents the single factor HJM 
model with stochastic volatility. The model is Markovianised using subsidiary state variables, 
as usual, and the specific volatility function enables a 3D PDE formulation. Section H] and 
[5] describe the full numerical method. The numerical validation is depicted in section |6l 
After some additional remarks regarding the numerical methods (section [7]), in section [S] we 
conclude. 

2. Interest rate market: the zero coupon bond 

In the interest rate market, the zero coupon bonds are taken as the most basic market 
instruments in the sense that any other financial quantity related to interest rates can be 
derived from them. These assets merely pay a monetary unit in a given future time that 
is called the expiration date, and the current price as a function of the expiration date is 
determined by the so called zero coupon bond curve (ZCBC). We will formally denote this 
curve as it follows: 

p{t,T),teR-^,T>t:T^p{t,-), (1) 

where t stands for the present time (valuation date), T stands for the bond expiration date 
and p (t, T) is merely the zero coupon bond price (ZCBP). We will also assume from now on 
that ZCBC fulfills the necessary regularity conditions. 
By definition, ZCBC is such that: 

• the ZCBP that expires at the present time is 1 (trivially): p{T,T) = 1. 

• ZCBP G (0, 1] and is a monotonically decreasing function (this property is assumed in 
order to avoid the existence of arbitrage). 

It is worth saying that in exceptional macroeconomic situations, it is actually possible that 
the latter property doesn't hold, due to the intervention of central banks: of Japan, in order 
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to encourage investors to buy their currency 21 



Provided the preceding properties, there always exists a function f(t,s),s > t such that 

,T>t 



f{t,s) > 0,Vs > t 



and consequently, / fulfills: 



f{t,T) 



dlnp{t,T) 
df 



VT > t. 



(2) 



(3) 



/^From a financial point of view, f{t,T) is interpreted in terms of the interest rate that an 
investor would receive if he sells in t a zero coupon bond with expiration date T and buys 
another with expiration date T + dT. This bond is the so called forward rate of interest. In 
particular, the short term rate of interest r fulfills: 



rit) = fit,t) 



d\np{t,T) 



dt 



(4) 



T=t 



and is of special interest because it quotes the return of an investment of one monetary unit 
in the present time, t, that is redeemed an infinitesimal time later, t + dt. 



3. HJM framework 

3.1. The model 

The single factor HJM model under hands is originally represented by the stochastic 
differential equation for the short rate r (details can be found in 0, Isl) 

dr{t) = (^^^ - <r{t) - /(O, t)) + y{t)^ dt + rjit, r{t))dW{t). 

dy{t) = [r]{t,x{t) f -2t^y{t)]dt, (5) 

where y{t) is a subsidiary state variable with no financial meaning, employed to Markovianise 
the model, k is a positive constant and dW {t) is a Wiener process. The volatility function 
rj{t,r{t)) is defined through 

r/(t,r(t)) = v^A(t)r(t)^W 

dv{t) = 9{1 - v{t))dt + e{t)^^)dZ{t) 

dZ{t) ■ dW{t) = pdt, (6) 

where dZ{t) is a Wiener process, X{t) and 7(t) are deterministic functions of time (7(t) G 
(0, 1]), v{t) is a stochastic variable that drives the rate variance, ^ is a constant that estimates 
the mean reversion speed of the process v{t), p is a constant that estimates the correlation 
between the short rate and the volatility (hence p G [—1, 1]), and e(t) is the volatility asso- 
ciated to v{t), which in this case is simply a deterministic function of time. 
Defining x{t) = r(t) — f{0,t), we have 

dx{t) = [-Kx{t) + y(t)]dt + r](t, x(t))dW{t), (7) 
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where dW (t) is another Wiener process. It can be analytically shown that the price of the 
zero coupon bond in this formulation reads 

G{s) = (8) 

K 

The HJM model is thereby fully characterized. The state variables belong to the following 
range: 

r(t) G [0,oo), v{t) e [0,oo), y{t) G [0, oo) 

3. 2. Model reformulation in terms of Partial Differential Equations 

Given the above model formulation, now we would like to price contracts whose future 
payoffs depend on the evolution of the yield curve, that is to say, those payoffs are deter- 
ministic functions of the ZCBP at certain times. Let C{t) be the price in a given time t, of 
a contract that pays in T > t in terms of F{T,p{T, s)), s > T, where P denotes the payoff 
function. The simplest example of this situation is the zero coupon bond whose payoff func- 
tion is a constant function, F{T,p{T, s)) = 1. The so called Caplet [H is another example of 
a payoff function, frequently traded in the market: 

F(T, p(T, Tm)) =max(l -Amp(T,Tm),0), (9) 

where T is the contract's expiration date, Tm is the contract's payment date and Am = 
1 + {Tm — T)K, where K stands for the strike and is a positive constant. In this paper we 
focus on both the Zero Coupon Bond and the Caplet as the derivatives under study. 
In these terms, the price of an interest rate derivative depends on time and space C{t, x, y, v) 
and fulfills the following partial differential equation: 

dC , d'^C , d'^C , d'^C dC dC dC 
ot or'' ov' orov or ov oy 



where the Feynman-Kac formula [12| is applied to map the stochastic problem into a PDE 
one (note that the financial notation for the derivatives -the so called Greeks- is not used 
here). The coefficients (rr, Cw, Cw, fJ^v, fJ^, fJ^y are functions of the space variables r,v,y and 
time, but their dependence has been omitted for clarity. These coefficients are given by 

Crr ^ lx{trr'<(% 
Cvv = ^e(t)^t; 

= A(t)r^We(t)pt; 
fir^^^^-<r-f{0,t)) + y 

Hy = 6{1 - V) 

= Xitfr^^^'^^v - 2Ky (11) 
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The boundary/terminal conditions of the problem are different depending on the payoff 
function. For the cases of the zero coupon bond and the Caplet, these are given by: 



• Zero coupon bond with expiration date T: The terminal condition simply reads 

C{T,r,y,v) = l (12) 
In the limit r — > oo the payoff is zero, and consequently one boundary condition is: 

C(t, r ^ oo, I/, = 0, t<T (13) 
Now, when r = the PDE reduces to 

CiT,r^O,y,v)^l, 
In the boundary y — > oo the PDE reads 

dC dC ^ dC ^ 

— - 2ny— = (15) 

C{T,r,y oo,v) = l 

that is analytically solvable (Lagrange's method yields a solution of the type C = 
F{ye'^^^,2r — y), where F is a generic function that must be determined to fulfill the 
boundary conditions). Finally, note that the boundaries v ^ oo and v — )■ are not 
relevant in this case as long as the zero coupon bond price is not a function of v. 

• Caplet: The terminal condition is given by 

C{T, r, y, v) = max(l - Amp{T, Tm] r, y, v), 0) (16) 
Just as in the case of the zero coupon PDE, the boundary condition for r — > oo is 

C{t, r ^ oo,y,v) = 0, t <T. (17) 

When r = the PDE reduces to: 

dC ^ . d^C dC . dC . dC , , 

^ + C..I.=o -^:;p^^'r^ + /^^l-o ^ + /^.l.=o - o (is) 

C(r, 0, y, v) = max(l - AMp(r, Tm; 0, y, v),Q) 

In the limit y ^ oo, the PDE and its boundary condition have the following shape: 

dC dC dC 

C(T, r, y — )■ oo, w) = 
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Thereby we have: 



C{t,r,y ^ oo,v)=0,t< T. 



(20) 



In the case w ^ oo we will have 



C{T, r,y,v ^oo)^ p{t, T) - p{t, Tm), 



(21) 



and when v — the PDE reduces to 




(22) 



C{T, r,y,v^O)^ max(l - Amp{T, Tm] r, y, 0), 0) 



that in this case can only be solved numerically. 

Additionally, when v — the following identity is commonly assumed to hold: 



and will be taken into account in the numerical development. 

3.3. Model parameters 

In order to fully specify the preceding PDE, the constants and functions implicitly defined 
in the HJM model have to be initialized. As a reference guide, and for the sake of an order 
of magnitude estimation, we introduce the characteristic market parameters in 2007: 

-Initial zero coupon curve: 



• K = 0.001 

• \{t) = 0.15 

• 7(i) = 0.9 

• e{t) = 1.5 

• e = 0.25 

. p = -0.75 



Observe that in a practical situation, for the calibration, values of the parameters need 
to be optimized using e.g. a conjugated gradient or genetic algorithm, iterating several 
times the model evaluation and comparing the results with the market data. This is also 
supposed to be done in real time, and therefore it is of fundamental importance that the 
model evaluation (numerical method) is as efficient as possible. In this work we focus on this 
fundamental issue and propose a numerical methodology that enables an efficient evaluation 
of the model (suitable for real time execution), while the global calibration problem (the 
aforesaid optimization method) is not addressed. 



• p{0,T) = 1.04-^ 
- Volatility function: 
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4. Optimizing the numerical scheme 

Before solving numerically equation ( !T0|) and in order to optimize the numerical approach 
and reduce its computational cost, we need to make some preliminary analysis. Three points 
are of major importance, namely (i) the study of some analytical solutions will provide in- 
formation about the system's solution itself, (ii) a detailed scale analysis of the problem will 
help to optimize the mesh resolution, increasing it only where/when needed. Finally, some 
considerations regarding the metric will also be addressed. 



4.1. Particular solutions 

The first case considers the solution of the zero coupon bond curve with expiration date 
T. As far as this solution does not depend explicitly on v (non-arbitrage conditions yield 
a volatility independent ZCBP), we assume v = without lack of generality in order to 
simplify the system of equations. Furthermore, constant k is generally small (volatility 
function depends on k) and we have thus assumed that it is also null. Hence, the partial 
differential equation reduces to: 

dC fdf{0,t) \ dC 



with the additional condition C{T,r,y) = 1. We can trivially map this PDE into a system 
of ordinary differential equations of the following shape: 

^^1, f^rC. (24) 

dT dr OT dr 

with initial conditions r = 0, t = T, r = ri, C = 1, and y being a parameter. 
It is indeed easy to check that its solution is: 

C{t,r,y) = e-/.^/(o,^)rf^e-(^-*)('-^(°'*»-^^. (25) 

Note that not assuming a null value for k is equivalent to substitute T — t by the function 
GiT — t), which is defined in (IHl). 

The usefulness of this analytical solution is twofold: first, it will serve to validate the nu- 
merical method, and second, it will stand as as boundary condition for the (more general) 
Caplet problem. 

4-2. Scale analysis 

In the financial realm, reducing the computing time as well as the computational cost 
(in terms of memory resource, for instance) is fundamental. An adequate temporal and 
spatial scale analysis will enable us to increase the mesh resolution only where/when needed, 
what leads to a saving of computational resources. Typically, this scale analysis is done by 
adimensionalizing the equations under study and consequently comparing the relevance of 
the respective terms (this technique is broadly used in fluid mechanics when performing the 



scale analysis of the Navier-Stokes equations, for instance 22|). 



Scales are indeed determined by the variables characteristic values, as well as by the boundary 
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conditions. In the case of the Cap problem, the variable of reference is the interest r(t) 
(r G [0, 100], that is to say, a percentage). We can rescale this variable as: 

f = r/ro, ro — 10^^ 

in such a way that the characteristic value of f is the unity (ro is usually the forward rate 
of interest observed at value date). Now, having in mind that the Cap's boundary condition 
reads 

C(T, 0, y, v) = max(l - (1 + (Tm - T)r^k)p{T, Tm] 0, y, t;), 0), 
where K = K/tq, we may define the characteristic time t as: 

i = rot. 

Finally, taking into account the relation between r and y (equation [24|) : 

Ar ~ yAt,-^y = y/r^. 
With these rescaled variables, the original equation takes the following shape: 

dC dC dC 

+h4{t,r,v,y)— + ^^5{v)-^ + he{t,r,v,y)— = rr^C (26) 

where 

hi{t,T,v) = ^A(t)V^Wro'^^(*^-^V;, h,{t,v) = ]^e{tfv, 

hiit, r, V, y) = (ro — ^ ^(r - /(O, t)) + roy), 

h{v)=e{l-v), h{t,r,v,y) = {\{tfr^^<^'\-2Ky). 

(Note that the "marks have been eliminated for notation simplicity). No rescaling has been 
applied to both v and C{t,r,v,y) as long as there's no dominant scale defined (the latter 
rescaling wouldn't affect the resultant equation). 

4.3. Metrics 

By introducing metrics in the independent variables (r, v, y), we can transform the prob- 
lem's domain into a computational domain which is usually simpler, and consequently con- 
centrate the mesh points in the areas under study. In our case of study (Cap), provided that 
the solution is likely to live in the strike's neighborhood (K), this should be the most dense 



zone. Following Tavella IJ], we have used an hyperbolic-like metric generally defined as: 

z = K + a sinh(c2X + Ci(l — x)). 
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ci = asinh ( — , C2 = asinh f , 

\ a J \ a J 

X e [0,1], z e [zq, Zc^], being z = (r, v, or y). (27) 

In order to get a mesh that is accurate enough in the zone under study, the parameters 
a, K, zq, z^ have to be correctly chosen for each dimension of the problem. 
It is worth saying that the introduction of these kind of metrics doesn't modify the system 
of equations in a substantial manner. Note that the derivatives with respect of the new 
variables can be expressed as: 

dC _dC 1 

dr dx dr/dx 

and 

^ ~ 'd^Jd^^Jdxf " 'diidrjdxf 

Thereby, if we denote Jz = dz/dx, and J2z = dz'^/dx'^, where 2; is a generic variable from 
(r, V, y) and x its respective transform, the original system of equations (|26|) will only differ 
from the new one in the substitution of functions hi, i = 1..5 by: 

gi{t,Xr,Xy) = hi{t,Xr,Xy)/J^, g2{t,Xy) = h2{t,X^)/jl, 

93 (^^-^r?-^!;) h^{t, Xj- , Xy) I J'l- j J^, Cj^ (t, X^ , Xy , Xy) /i4(t, Xj-, X'q, Xy^ j J^- ll-i{t, X^ , X^) J^r I J x "i 
y^it, Xy^ hc^{Xy) I J J, h2(t, Xy)J2v/ J^]) dsity -^ri -^v: -^y) hgit, Xj-, Xy, Xy) j Jy. 

While hitherto we have only defined metrics in the direction of the independent coordinates, 
it is actually possible to employ more complex transformations that involve several variables 
= ^{fyV,y)). However, it is likely that the integration domain wouldn't in that case be 
cartesian anymore, and consequently the finite difference scheme wouldn't apply. Further- 
more, while these types of metrics would eventually enable us to eliminate the cross derivative 
terms in (12^ (transforming the original equation into its canonical form), this transforma- 
tion would on the other side modify the frontiers of the problem from straight to curve lines, 
something that not desirable in any case. 

Once the preliminary insights have been put forward, we will describe in the next sections 
the numerical method employed to integrate equation f l26l) as well as the results that we have 
obtained. 



5. Numerical methods 

As commented above, the equation under hands is (!26|) : 

dC , ,d^C , , d^C , , d'^C 
'o^+9iit,r,v)^ + g2it,v)^ + gs{t,r,v)^ + 

dC dC dC 

+94{t,r,v,y)— + g5{v)— + gQ{t,r,v,y)— = ttqC. {2i 
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Let us define a new temporal variable 

t = T -t 

where T stands for the maturity. The integration should then be done for t G [0,T]. Note 
that equation ( 126|) is parabolic for r and v and hyperbolic for y. 

In the particular problems concerning the estimation of financial derivatives, the execution 
time of numerical tools is an issue of fundamental importance. According to this fact, it seems 
suitable to apply second order schemes for the discretization of both temporal and spatial 
partial derivatives, as far as these schemes show optimal computational cost and adequate 
precision. In a second step, one has to decide whether to apply explicit or implicit schemes. 
Naturally, the simplest option is always to tackle explicit schemes, which are fast and easy 
to implement. However, it is easy to check that due to the second derivatives, the following 
relation holds for the temporal and spatial resolutions 

At ~ {Arf. 

This relation implies that in order to achieve a precision of say 10~^ in the solution, a 
time step would need 10"^ iterations. Moreover, the coefficients of the derivatives are powers 
of r, V or y, and the integration domain ranges to the infinite. Since the time step is inversely 
proportional to those coefficients, the problem comes to be even more delicate. We can thus 
conclude that explicit temporal schemes won't fit in this case due to their inevitably lengthy 
behavior. Thereby, we will have to choose implicit schemes for the temporal integration. 
These have the following general expression: 

dC 

— = F{C) -> - = At (0F(C"+i) + (1 - ^)F(C")), 

where we have employed the usual numerical methods notation C" = and t" = nAt. 

Here 6 stands for an explicit scheme for = (Euler scheme) while it stands for implicit 
schemes when 6^0. More concretely, 6 = 1 characterizes the so called Euler implicit scheme 
and finally 6 = 1/2 characterizes the second order Crank-Nicholson scheme. We will use the 
latter one as the temporal integrator as it is adequate to be used within ADI schemes (this 
will be explained further in the text). 

The first, second and cross spatial derivatives, are discretized by centered finite difference 
schemes as it follows: 



^Xx'^ r\ o" A o 1 ^ X^^ r\ T 7 

ax"' Ax^ ox Ax 

OxzU = o o = ^ — 7. — (29) 

oxoz AxAz 

where x, z is a generic variable that stands for r, v or y. The notation 5xxU., SxU and SxzU 
describes the above difference schemes. Gathering both spatial and temporal schemes, we 
come to a final discretization of the following kind: 

^n+i _jjn ^ At(^F(f/"+i) + (1 - 6)F{U'')) (30) 
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where 



F{U) = gi{t,r,v)6rrU + g2{t,v)6^M + g3{t,r,v)6rvU 
+gi{t, r, V, y)5rU + g^{v)5^U + g^^it, r, v, y)SyU - rroU, 



and U = {Cijk,i = 0..nr,j = 0..nv,k = 0..ny}, (31) 

is the discretized solution vector in a structured mesh of dimension {nr,nv,ny). 

It is worth saying that the use of centered difference schemes allows us to obtain a compact 

stencil. For instance, note that the discretized equation in the point Cijk only contains 

information of {Ci+ij^k, Cijk, Ci^ijk} ■ Focusing on variable r, equation fl5U]) would adopt the 

shape: 

jjn+l _Jjn ^ At{{eLrU"+^ + 6) + (1 - 0)(L,f/" + 6)), 



where Lr is a tridiagonal matrix [13| representing the spatial discretization of gi{t,r,v)6rr + 
g4{t,r,v,y)6r in F{U) (equation (1ST]) ), according to the discretization schemes depicted in 
equation fl29l) . L^, Ly and L^v will be defined equivalently (see below). 
Tridiagonal systems are indeed quite easy to implement and solve (for instance, the Thomas 



algorithm [23[ solves a tridiagonal system in 7A^ operations, where N is the order of the 
system) . This goodness will enable the use of Alternating Direction Implicit schemes (ADI) 
YlA isl as will be shown further in the text. 



Finally, taking into account that the equation is indeed linear, it can be written as: 

f/"+i -U- = 9 At [L,(r,r„i;„yfc)t/"+i + L,{e,v,)U^^^ + L,(r, r„ y,)f/"+i] 
+ (1 - 9) At [L,(r, r„ V,, yuW + ^.(^", v,)U'' + L,(r, r„ v,, y^W] 

+AtL™(r, Ti, Vj)U'' + AtBie, Ti, Vj), 

and realigning, 

[/ - 9At{Lr{t''+^/\ r„ v,,yk) + L,{t^^^/\ v,) + Ly{e+^/\ r„ v,, yu)] (f/"+' - U^) 
= At [L,(r, Vj, yk) + L„(t", Vj) + Ly{r, Vj, yu) + L™(r, v^)] f/" + At5(r, r^, Vj) 

= AtF(?7",r,r„i;„yfe).(32) 

Note that the operators Lr,Ly,Ly include the terms related to the spatial discretization. 
These operators, treated implicitly, give rise to a system of equations(/ — 6'(Lr + L^ + Lj^))f/ = 
/) which in general has 7 diagonals, and whose resolution can be performed applying either 
direct or iterative methods. However, each one of them treated separately can be rewritten 
as a tridiagonal matrix, whose resolution is trivial as commented above. Note also that the 
operators Lr,Ly,Ly have Neumann boundary conditions and consequently do not include 
any Dirichlet-like information. 



Finally, note that the mixed derivative term (L^^,) is only treated explicitly, because 
otherwise its inclusion in the implicit scheme would eliminate the tridiagonal structures, and 
would consequently avoid the use of ADI schemes that will be described in what follows. This 
fact does not affect in any case neither to the convergence nor the precision of the numerical 
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solution. Detailed numerical analysis and validation of the mixed derivative term has been 
already performed by different authors: examples of implementation for the fluid mechanics 



convergence can be found in 



problems are given in [25|, |26|. and a detailed analysis of numerical stability analysis and 
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5.1. Alternating Direction Implicit (ADI) schemes 



29 , used in the 



The ADI schemes belong to the category of Splitting methods |2J, |28|, 1 131 . 
resolution of multidimensional PDE systems. The key idea behind these methods is to sepa- 
rate the original multidimensional problem in several unidimensional split problems. Then, 
each split problem can be under certain conditions reduced to the resolution of a tridiagonal 
system of equations. These conditions are related to the use of centered spatial operators in 
structured meshes, which is our case. 

in order 
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ADI schemes were initially introduced by Douglas, Peaceman and Rachford |24J . 
to integrate, using finite difference schemes, the well known Navier-Stokes equations describ- 
ing the fluid motion. Some modifications have been put forward so far (see for instance 
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25l |27| . Il3l | , in order to apply these schemes to either stationary or non stationary prob- 
In this work we will use an ADI scheme recently put forward by Hout & Welfert liT 



lems 

called the Douglas scheme 



Consider equation f l5^ . this one can be formally written as: 

[/ - eAt{Lr + L^ + Ly)]AU'' = AtF(f/") 
The Douglas scheme applies thus in the following way: 

Af/° = AtF(f/") 
[I - At Lr AU^] = Af/° 
[/ - At Lv AU'^] = AU^ 
[I - At Ly Af/"] = Af/2 

jjn+l ^ jjn ^ ^fjn 

Note that each step only requires the resolution of a tridiagonal system of dimension nr, 
nv, or ny. The unconditional convergence of this scheme has been proved for 9 = 0.5 (Crank- 
Nicholson) in 2-dimensional systems with constant coefficients jl3|. However no similar study 
has been performed so far in the 3-dimensional case with variable coefficients i^l, which is 
nonetheless our case. Special attention will be thus paid to the convergence behavior of the 
solution. 



It is worth saying at this point that the Craig & Sneyd [30[ method is an apparent 
improvement to the Douglas scheme (in terms of the solution precision) when mixed derivative 
are present, while being more expensive computationally speaking. We actually have also 
tackled this ADI scheme, but given that no such improvement has been observed, we will 
only focus on the Douglas scheme. 

Finally, the Douglas scheme, as any other ADI scheme, is an Approximate Factorization 
(AF) of the original equations with the errors of order At^ for 3D problems. An efficient 
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subiteration procedure can be applied to eliminate the AF. This method, known as Huang's 
approximate factorization correction 
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has been also checked in this context, but it is 
computationally more expensive and no additional improvements have been observed in the 
range of accuracy we are working. 



6. Validation 

We have done two kind of studies in order to validate the numerical methods: 

(i) first, we have compared the numerical solution of the zero coupon curve or deduction 
curve with its analytical solution, in two different situations, and 

(ii) second, we have compared the numerical solution of a Caplet with the one obtained by a 



2D-Heston model 13j (it is easy to check that the model under study behaves, for Tm — )■ T, 



as a Heston model for the libor rate with identical parameters). 

6.1. Zero coupon curves 

In the first study we consider the parameters and the initial zero coupon curve depicted 
in section [3731 Thereby, the forward rate of interest is constant f{0,t) = log(1.04) and it has 
a null derivative. In figures we compare the theoretical zero coupon curve with the one 
obtained through numerical simulations with r = log(1.04), y = 0, and t = T. Concretely, 
figure ([1]) shows the error for different meshes. Note that this one is always below 10^^ even 
for coarse meshes. As a result, we have set the mesh reference values to 100 x 40. In figure 
([2]) we plot the convergence of the solution as a function of the number of time steps per year. 
Notice that from 12 steps per year, in a given mesh the variations are quite small (O(10~^)). 



In figures ([3]) and (j4]) we show the error's spatial distribution for a given set of parameters, 
assuming yoo = 25 (figure ([3])) or ?/oo = 250 (figure (jlj)) respectively. Note that in the 
former case, some non desirable errors take place in the infinite boundary, which can actually 
propagate into the zone under study y ~ (figure ([2D). 

Up to know the numerical method is validated, as far as the solution's error is confined, 
in the zone under study, around 10^^. However, as long as the analytical solution strongly 
depends on the initial curve, it is necessary to check whether if the precision of the numerical 
method holds for more realistic curves (with non null forward rate interest curve derivative). 
For that task, in a second example we tackle a new initial zero coupon curve, which is not 
anymore a continuous curve but a discrete valued one (figure (jS]))- Its derivative is plotted 
in figure ([6]) and stands for the forward rate of interest curve, and its second derivative is 
plotted in figure ([7]). 



As long as the curve is expressed in terms of discrete values, we need to perform a 
smoothing approximation in order to introduce it in the simulation. Notice that the solution's 
smoothness will strongly depend with the smoothness of this initial curve (this is due to the 
fact that the temporal derivative of the solution is related to the derivative of the forward 
rate of interest curve). An appropriate solution to this problem is to approximate the initial 



curve with splines [23[, in order to have a piecewise function with continuous second derivative 



df{0,T)/dT, and consequently have a solution with continuous temporal derivative. 
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According to this approximation, we have performed the same simulations and analysis 
as for the first study. Conclusions are plotted in figures (IHlfTOl). As expected, the fact that 
the derivative of the forward rate of interest curve is non null has a net effect in the precision 
of the solution. While it is quite easy to achieve convergence of order O(10~^), it comes 
necessary to overrefine the mesh (figure ([H])) or alternatively increase the number of time 
steps (figure (^) in order to go beyond 10~^. Nevertheless, as is shown in figure (ITUI) . the 
error's spatial distribution is quite similar to the one found in the first study: we can conclude 
that the numerical method correctly reproduces the expected results. 

6.2. Caplet 

The second validation test consists in making a comparison between the results obtained 
with several Caplets and those obtained by the Heston model, which is an already validated 



model 13 



In order to optimize the mesh's size, we have performed a previous analysis of the numerical 
scheme's convergence (both in spatial and temporal discretizations). Some of the results 
are plotted in figure ( TTT]) . where we represent the evolution of a generic Caplet's prime as a 
function of nr, nv, ny and nt. Notice that we need at least a mesh size of 100x40x40 if we 
seek variations of the prime below 10~^. With a reference mesh of 100x40x40, the number of 
time steps does not affect practically the results. 

The difference between the solutions of the two models are plotted in figures f[T^ and 
f[T^ . both for the premium and for the volatility. 

6.3. Computation time 

In the following table we have plotted, as a reference guide, the required computation 
time for different Caplets. Simulations have been run in a mesh of 100 x 50 x 50, with 
nt = 12 steps per year, in a Pentium{R)IV processor (3.2 GHz, 1Gb RAM). Results are 
quite satisfactory. 



TMc 


Tc 


Ntotal 


CpuTime 


2 


1 


12 


0.9 s 


11 


10 


120 


7.3 s 


20 


19 


228 


14.0 s 



For illustration, the numerical solutions obtained for the Caplet's prime and the greeks 
P = ^ and vega, = ^ are shown in figures fll7lll8|19p respectively. The computations 
have been performed for TMc=2, Tc=l and Ntotal=12 and the pictures are shown at y=0. 



7. Some additional numerical aspects 
7.1. Metric choice 

As commented in section 14.31 it is highly recommendable to introduce a metric layer in 
the numerical method, such that the domain under study transforms into a computational 
domain which is typically easier to handle (in most cases, this one is the unity cube), as long 
as this domain enables the use of structures uniform meshes, where one can concentrate the 
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mesh points wherever needed. The mesh that has been used in this work is hyperbohc (see 



eq. [27|) . following Tavella One of the main properties of these meshes is that one can 
concentrate as many points as needed in the inner regions of the zone under study, in order to 
achieve a better resolution. It is thus convenient to fix the parameters related to the domain 
transformation. For instance, r will transform according to: 

r = Kr + ar sinh(c2r-a;r. + Cir(l — Xr)), 

'{ro-Kr)\ . f{roc-Kr] 

Cir = asmh , C2r = asmh 



Otr J \ ar 

where the jacobian of the transformation reads: 

-J— = ar Smh{c2Xr + Ci(l - Xr)){c2r - Clr). 

We have then four parameters Kr,ar,rQ, to fix: 

• '"o, ''^oo define the real domain of study. Obviously ro = 0. On the other side, r^o must be 
such that his values doesn't modify the solution in the zone under study (that is, close 
to the strike). There is no recipe in order to find the adequate value, but after some 
preliminary estimations and taking into account the boundary conditions, we have set 
Too = 250. 

• Kr defines the region under study, that is, a neighborhood of the strike. 

• ar This parameter provides a measure of the mesh's stretching, i.e. the number of 
points that will be concentrated in the zone of interest -close to the strike-. Concretely, 
the smaller a^, the larger concentration. Given that its value also affects the jacobian 
of the transformation, it is desirable that ar is such that the jacobian be close to 1. In 
figure (fT4l) we plot this dependence, for r = l(ro). Note that for ~ 0.05 — 0.1, the 
jacobian reaches the unity. A similar study forvey lead us to fix ay ~ 0.5 y ~ 0.05. 

As a summary, the metric's characteristic values are: 

ro = Vo = yo = 0, = = 250, v^o = 30, 



ar = ay = 0.05, Kr = Strike, Ky = 0, = input{c^ 0.5) 
7.2. Softening of the initial condition 



Following Tavella [IJ], as far as the payoff is typically a discontinuous function, small 
variations in the strike lead to a non smooth behavior of the solution. This is not desirable 
and therefore some numerical techniques should be applied in order to soften it: 

• Perform a dynamical modification of the mesh, related to the payoff's shape. This is 
an elegant solution, however for practical purposes this technique is not well fitted as 
long as it usually leaves to mesh interpolation. 
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• Soften the initial/final conditions. In order to do so, one can define an average ini- 
tial/final condition in the following terms: 

Cijk{t,r,v,y) = - [ C{t,r,v,y)du, 

where a; is a control element centered in the mesh point ri,Vj,yk- The effect of this 
average is represented in figures fllStllGp . 

7. 3. Boundary condition for variable y 

Note that the HJM model under hands is only convective for variable y (first derivatives 
are null for every variable but y). When y = 0, the solution's characteristic crosses the 
domain, what indicates that the boundary takes some information from inside. Consequently, 
the discretization of both the interior and the boundary should be consistent, and then a 
second order scheme should be applied to the boundary discretization. We have implemented 
two different possibilities in the numerical scheme: 

• Advanced first order differences: 

dy ~ Ay 

• Advanced second order differences: 

dy ~ 2Ay 

and quite surprisingly, no significative differences have been found between both schemes 
results. 

8. Conclusions 

In this paper we have proposed a complete numerical methodology to efficiently solve a 
single factor HJM model with stochastic volatility. For this task we have first Markovianised 
and reformulated the model in terms of a three dimensional PDE system. The numerical 
method involved finite-difference and Crank-Nicholson schemes for the spatial and temporal 
discretization respectively. In order to decrease the computing time without loosing preci- 
sion, we have successfully applied ADI schemes and performed a preliminary scale analysis 
to optimize the mesh resolution. The validation of the numerical schemes has been done 
comparing the numerical solution of some test curves (zero coupon bond, Caplet) with ana- 
lytical models and a Heston model respectively. The goodness of the results in terms of low 
computation time (order of seconds in a standard pc) and good accuracy (typical errors ob- 
tained either for artificial or quite realistic forward interest rate curves haven't gone beyond 
10""^ in any case) suggest that the method is suitable to be applied in realistic applications, 
concretely in the financial industry. 

Acknowledgments The authors thank anonymous referees for their helpful suggestions. LL 
acknowledges financial support from grants FIS2009-13690 and S2009ESP-1691. 



17 



References 

[1] John C. Hull, Options, Futures, and Other Derivatives, 4th Edition, Prentice Hall, 
(1999). 

[2] R. Rebonato, Interest-Rate Term-Structure Pricing Models: A Review, Proc. Roy. Soc. 
Lond. A 460, 2043 (2003). 

[3] Heath, Jarrow, Morton, Bond pricing and the Term structure of interest rates: a new 
methodology for contingent claim valuation, Econometrica 60, 1 (1992). 

[4] L. Andersen and V. Piterbarg, Interest Rate Modeling, Atlantic Financial Press, (2010). 

[5] P. Ritchken and L. Sankarasubramanian, On Markovian Representations of the Term 
Structure, Working Paper series 9214, Federal Reserve Bank of Cleveland (1992). 

[6] M.W. Baxter, General interest rate models and the Universality of HJM, in Financial 
Derivatives, edited by M.A.H. Dempster and S.R. Phska, CUP (1997). 

[7] A. Caverhill, When is the short rate markovian?. Mathematical Finance 4, 4 (1994). 

[8] P. Ritchken, L. Sankarasubramanian, Volatility structures of forward rates and the dy- 
namics of the term structure. Mathematical Finance 5, 1 (1995). 

[9] K. Inui, M. Kijima, A Markovian framework in multi-factor Heath-Jarrow-Morton mod- 
els. Journal of Financial and Quantitative Analysis 33, 3 (1998). 

[10] C. Chiarella, O.K. Kwon, Forward rate dependent Markovian transformations of the 
Heath-Jarrow-Morton term structure model. Finance and Stochastics 5, 2 (2001). 

[11] C. Chiarella, O.K. Kwon, Finite dimensional affine realisations of HJM models in terms 
of forward rates and yields. Review of Derivatives Research 6, 129-155 (2003). 

[12] I. Karatzas, S.E. Shreve, Brownian motion and stochastic calculus. Graduate texts in 
Mathematics, 2ed, Springer (1991). 

[13] D.J. Duffy, Finite Difference methods in financial engineering: a partial differential 
equation approach, Wiley (2006). 

[14] D. Tavella, C. Randall, Pricing Financial Instruments: The Finite Difference Method, 
Wiley (2000). 

[15] O. Cheyette, Markov Representation of the Heath-Jarrow-Morton Model. Available at 
SSRN: http://ssrn.com/abstra ct=6073| or doi:10.2139/ssrn.6073, March 26, Barclays - 
San Francisco, Ca Office (2001). 

[16] Paper: V. Piterbarg , TARNs: Models, Valuation, Risk Sensitivities, Wilmott magazine 
(2004) pp: 62-71. 



18 



[17] R. Rebonato, On the simultaneous calibration of multi-factor log-normal interest-rate 
models to Black volatilities and to the correlation matrix, Journal of computational 
finance 2, 4 (1999). 

[18] R. Rebonato, Volatility and Correlation, 2nd edition, John Wiley & Sons (2004). 

[19] P. Cotton, J-R Fouque, G. Papanicolaou, R. Sircar, Stochastic Volatility Corrections 
for interest rate derivatives. Mathematical Finance 14, 2 (2004). 

[20] C. Chiarella, O.K. Kwon, A complete Markovian Stochastic Volatility model in the HJM 
framework, Asia-Pacific Financial Markets 7, 4 (2000). 

[2 1] http: / / money . cnn.com /1998/11/06/ economy / j apanhank / . 

[22] Munson, Fundamentals Of Fluid Mechanics, Wiley (2007). 

[23] The book Numerical Recipes: The Art of Scientific Computing, Cambridge University 
Press (2007). 

[24] D.W. Peaceman, H.H. Rachford, The numerical solution of parabolic and elliptic differ- 
ential equations, SIAM 3, (1955) pp. 28-41. 

[25] S. McKee, D.P. Wall, and S.K. Wilson, An Alternating Direction Implicit Scheme for 
Parabolic Equations with Mixed Derivative and Convective Terms, Journal of Comput. 
Physics. 126, 64-76 (1996). 

[26] G. H. Klopfcr, R.F. Van dcr Wijngaart, C.M. Hung and J.T. Onufcr, A Diagonalizcd 
Dominant Alternating Direction Implicit (D3ADI) Scheme and Subiteration Correction, 
AIAA-1998-2824 (1998) 

[27] K.J. Hout and B.D. Welfert, Stability of ADI Schemes Applied to Convection- Diffusion 
Equations with Mixed Derivative Terms., Applied Numerical Mathematics 57, 1,19-35 
(2007). 

[28] J. Douglas, H.H. Rachford, On the numerical solution of heat conduction problems in 
two and three space variables. Trans. Amer. Math. Soc. 82, (1956) pp. 421-439. 

[29] J.L. Steger, Implicit Finite Difference Simulation Flow About Arbitrary Geometries with 
Applications to Airfoils, AIAA Paper 77-665 (1977). 

[30] I.J.D. Craig and A.D. Sneyd, An Alternating-Direction Implicit Scheme for Parabolic 
Equations with Mixed Derivatives, Comp. Math. Appl. 16 (1988). 



9. Figures 



19 




Maturity [years] 

Figure 1: Plot of the numerical solution's error as a function of the mesh size, for the initial zero coupon 
curve. The particular values of the metric are — 0.05, ay — 0.5, Too = 25, yoo = 250. 
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Figure 2: Plot of the numerical solution's error (convergence) as a function of the number of time steps per 
year, for the reference mesh of size 100x40. The metric parameters are the same as for figure 1. 
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Figure 3: Error spatial distribution obtained in t=T=20 years for the zero coupon curve. The metric 
parameters are — 0.05, aj, — 0.5, Too = 25, 2/00 = 25., the mesh is the reference one and the temporal 
discretization assumes 12 time steps per year. 
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Figure 4: Error spatial distribution obtained in t=T=:20 years for the zero coupon curve. The metric 
parameters are ar — 0.05, — 0.5, Too — 25,?/oo — 250., the mesh is the reference one and the temporal 
discretization assumes 12 time steps per year. 
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Figure 5: Initial zero coupon curve (discrete curve). Note that the plot is in semilog. 
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Figure 6: Approximation of the forward rate of interest curve, obtained with splines. 
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Figure 7: Approximation of the forward rate of interest curve derivative, obtained with sphnes. 
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Figure 8: Plot of the numerical solution's error as a function of the mesh size Nr x Ny, for the initial zero 
coupon curve. The particular values of the metric are = 0.05, ay = 0.5, Too — 25, yoo — 250. 
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Figure 9: Plot of the numerical solution's error (convergence) as a function of the number of time steps per 
year, for the reference mesh of size 100x40. 




Figure 10: Error spatial distribution obtained in t=T=20 years for the zero coupon curve. The metric 
parameters are = 0.05, aj, = 0.5,7-00 = 150,?/oo = 150., the mesh is the reference one and the temporal 
discretization assumes 12 time steps per year. 
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Figure 12: Difference between the Heston model and the HJM model in the assessment (premium) of several 
Caplets. Metric parameters: = 0.5, = 0.5, — 0.5, Too = 250, Woo ~ 30,?/oo = 250. and 12 time steps 
per year. 
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Figure 13: Difference between the Heston model and the HJM model in the assessment (volatility) of several 
Caplets. Metric parameters: ar = 0.5, a« = 0.5, ay — 0.5, roo = 250, Uoo — 30,2/oo = 250. and 12 time steps 
per year. 
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strike 

Figure 15: Premium variation for a generic Caplet and different strikes. Simulations have been realized in a 
single mesh of 100x40x40 with nt=12 steps per year. 
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Figure 16: Variation of the premium derivative for a generic Caplet and different strikes. Simulations have 
been realized in a single mesh of 100x40x40 with nt=12 steps per year. 
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Figure 17: 3-Dimensional view of the premium for a generic Caplet at y=:0. Simulations have been realized 
in a single mesh of 100x50x50 with nt=12 steps per year. 
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Figure 18: 3-Diniensional view of p for a generic Caplet at y=0. Simulations have been realized in a single 
mesh of 100x50x50 with nt=12 steps per year. 




Figure 19: 3-Dimensional view of vega for a generic Caplet at y=0. Simulations have been reaUzed in a single 
mesh of 100x50x50 with nt=12 steps per year. 
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